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Aims. Collections of dust, grains, and planetesimals are often treated as a pressureless fluid. We study the validity of neglecting the 
pressure of such a fluid by computing it exactly for the case of particles settling in a disk. 

Methods. We solve a modified collisionless Boltzmann equation for the particles and compute the corresponding moments of the 
phase space distribution: density, momentum, and pressure. 

Results. We find that whenever the Stokes number, defined as the ratio of the gas drag timescale to the orbital timescale, is more than 
1 12, the particle fluid cannot be considered as pressureless. While we show it only in the simple case of particles settling in a laminar 
disk, this property is likely to remain true for most flows, including turbulent flows. 
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1. Introduction 

The interaction of solid particles with gas in protoplanetary 
disks is one of the key mechanisms for planetesimals and 
planet formation. Moreover, this fundamental process gets more 
and more attention for other astrophysical media (interstellar 
medium, planetary atmospheres, etc.) as grains can dampen 
fluid motions and be the targets of adsorption and condensation 
of molecules and can be the base for complex chemistry, e.g., 
the formation of the molecule H2 in the interstellar medium. 
With the advance of observational facilities and techniques, dust 
dynamics is understood more and more as a crucial process 
for most media, while, due to the discrete nature of grains, 
modelling their dynamics remains tricky. 

There are two main methods for computing the dynamics 
of particles embedded in a fluid. One can either study a 
large number of individual particles using N-body techniques 
or consider the particles as a fluid, which is treated with a 
modified Navier-Stokes/Euler equation. While the first method 
is reliable and precise, it can only be applied to a finite (and 
in fact limited) number of particles; therefore, macroscopic 
quantities reconstructed from the N-Body simulations outputs 
are sometimes inaccurate due to a low sampling of the phase 
space. The second method enables one to consider an almost 
infinite ensemble of particles by computing their macro- 
scopic properties (density, velocity field, etc.). This makes the 
"two-fluid" approach (gas and solids) the preferred method of 
treating the dynamics and growth of solid particles in a fluid, 
especially in the planet formation community. Consequently, 
a dense and rich literature is based on this approach, and 



many fundamental m echanisms of planetary formation have 
been unraveled (e.g. ICuzzi et ail 1 1 9931: iDubrulle et al.l Il995; 
IStepinski & Valageasl 119961 jBarriere-Fouchet et al. 2005 
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llnaba et al.l 120051: iFromang & Papaloizou 2006; Johanse n et al.l 
l2006h . 

The two-fluid approach is usually implemented by introduc- 
ing two Euler equations, coupled through the gas drag force. 
Moreover, in the fluid equation for particles, the pressure term 
is usually ignored. This approximation is often justified by the 
statement that collisions between particles and gas are much 
more probable than collisions between particles making the 
pressure term negligible compared to the gas drag force. 

The purpose of the present paper is to question the validity 
of ignoring the pressure of the particle fluid by computing 
it in the very simple case of particles settling in a disk. The 
vertical dynamics of a solid body in a disk can be understood 
as a damped oscillator. In a Keplerian disk, if gas is absent or 
very tenuous, the particle oscillates around the midplane at a 
frequency equal to the orbital frequency (the oscillation comes 
from inclined orbits). The friction with the gas damps this 
oscillation, as the particle feels a "headwind" when it crosses 
the gas disk. Like any damped oscillator, there are two regimes: 
an underdamped regime for low damping and an overdamped 
regime, where the particle has no time to oscillate for a full 
period. This occurs for high damping rates. 

We can compute the pressure of the solid fluid embedded 
in a laminar gas disk by solving exactly the collisionless 
Boltzmann (or Vlasov) equation for the phase space density of 
particles. As we show, as soon as the particles are larger than 
some critical size, the pressure cannot be neglected. This work 
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extends the work of iGaraud et alj d2004h . who comes to similar 
conclusions based on asymptotic calculations. 

In Sect. 2 is developed the general formalism. Section 3 is 
devoted to a very simple test case that can be used to check the 
validity of more general solutions and intuit how the pressure 
varies. In Sect. 4, the exact general solution of the Vlasov equa- 
tion is given. In Sect. 5, this solution is applied to a special case 
that retains most of the expected general properties. A surpris- 
ingly simple solution is found for the pressure and the scaling 
laws, and asymptotics are discussed. In Sect. 6, we discuss our 
results and their implications for more complex and turbulent 
flows and give our conclusions. 

2. General formalism 

We consider the equation for the evolution of the phase space 
density in 1 + 1 dimensions (one in space, one in velocity) under 
the action of a drag force from an inert fluid. This is not exactly 
the common collisionless Boltzmann or Vlasov equation, as the 
drag force is not conservative. Instead, the general Vlaso v equa- 
tion with non-conservative forces takes the form (see e. g.[Lutsko 
fl997tlCarbalhdo et alJl2006t lYoudin & LithwicU2007l) : 

M+f(-/H <» 

at oz ou \ m I 

where z is the space coordinate, u is the velocity coordinate, F/m 
the non-conservative force per mass (acceleration), and / the 
phase space density (or probability density function). When / 
is known, the macroscopic quantities (density, velocity, energy) 
are moments of the distribution /. In particular, the zeroth order 
moment is the density n(z, t): 

X+oo 
f(z, u, t)du, (2) 

the first-order moment is the momentum (the product of the den- 
sity by the macroscopic velocity 5) 

r»+oo 

n(z, t) u(z, t) — I uf(z,u,t)du, (3) 

\J — oo 

and the second-order moment is related to the pressure (a scalar 
in ID): 



path of the underlying gas. It should, however, be noted that 
the mean free path l m f p in a circumstellar disk is very large, 

typically varying like l m f p = 10 8 ( iooau ) ~ cm ' assumm g that 

the gas number density n g varies like n g = 10 8 ( 1(y f AII ) ' cm -3 
dGuilloteau & Dutrevlll998l:lHersant e7ai]|2005l) . 

Under these conditions, Eq. ([1) becomes 

d l f + ud z f-d u (co 2 z + ^ ) )f = o (5) 

where z is now the vertical space variable and u the vertical ve- 
locity. To simplify the notations, we use the convention d x — 4-. 
This equation can be developed into pure advective terms: 

8,f + udj - (a?z + dj - t = (6) 

The last term on the lefthand side (— f/r) is the difference 
with a Vlasov equation for conservative forces and deserves 
some comments. This term induces an exponential growth of 
/. Without solving the equation, one can understand that the 
conservation of the number of particles (following from Eq. [TJ, 
requires that the surface of the phase space field of density / has 
to collapse in one direction so that, at infinite times, the field of 
/ tends to a discontinuous distribution, consisting of a family of 
curves, but nothing like a continuous surface. This property of 
the equation remains valid even in 3+3 dimensions and with gas 
turbulence. 

At this point, it is convenient to adimension the time, intro- 
ducing two new variables t' = tot and «' = u/a>. Using these new 
variables and omitting the primes for clarity, we get 

d,f + ud-J-(z+^)d u f- f -=0 (7) 

where S = cot is a Stokes number of the problem. For small 
(large) particles, S is small (large) due to the strong (weak) 
coupling of the particles with the fluid. 

We now can artificially get rid of the exponential growth 
term by defining 



f + °° 2 

P(z, t) = I (u - u) f(z, u, i)du 

xj — OO 



* = fe 



-i/s 



f 



u 2 f(z, w, t)du — nu 2 . 



(4) 



In the case of particles settling in a disk, there are two 
forces: gravity and gas drag. When the disk is geometrically 
thin, the vertical gravity can be written as £ = -co 2 z, where 
to is the orbital frequency, while the drag force, in the absence 
of gas motions, can be written as £ = — u/r where r is the 
so-called stopping time dWhippld 1 1972b . This is the charac- 
teristic timescale for a particle to stop in a fluid. Here, we 
consider r as a free parameter, which would in general depend 
upon the particle size, morphology, density, the density of the 
gas, and the turbulent intensity in the grain trail. Here, we 
consider it as independent of both space and velocity. This 
is an important simplification, corresponding to the so-called 
Epstein regime, where grains are smaller than the mean free 



and finally get 



(8) 



(9) 



The behavior of this equation is rather easy to foresee. For in- 
finite S , the remaining terms correspond to a rotation in phase 
space. This rotation is a direct consequence of energy conserva- 
tion. In real space, particles oscillate around the midplane with 
zero velocity at the highest absolute altitude and maximum ve- 
locity when particles cross the midplane. The gas-drag term in- 
duces a collapse along the velocity axis. The net outcome of both 
the rotation and the collapse along the velocity axis (which cor- 
responds in the real space to the dampening by gas drag of the 
inclination of the particle orbit) is a global collapse towards the 
center of the phase space (z = 0, u = 0), but not a spherical 
collapse. 
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3. Simple case 



5. Particular case for the initial conditions 



Although the general solution is given in the next section, it is 
useful to solve Eq. ((9) in the simplest case: infinitely large or 
massive particles (S = +oo). In this case, the equation admits a 
simple stationary solution with separated variables: 



Computing the moments of the phase space distribution cannot 
be done without the initial conditions. Here, we choose an initial 
condition composed of Gaussians, with one Gaussian per dimen- 
sion, which is relevant to most physically meaningful cases: 



i// = Ke' 



(10) 



Using the definition of the moments, the pressure is straightfor- 
ward: 



P(z) = n(z)cr 2 . 



where the density is 



n(z) = Kcr ^l2ne 



(ID 



(12) 



Equation (fTTb is similar to an equation of state for isothermal 
fluids, where cr plays the role of a (constant) sound speed. As we 
will see in Sect. [5] this behavior remains valid even in more gen- 
eral cases. With the equation of state (fTTb . solution ( fT2l is the 
result of the hydrostatic equilibrium of an isothermal gas, the 
pressure force balancing the vertical gravity. The stationary so- 
lution (fTOl being even in u, it corresponds to a zero macroscopic 
velocity. In this case, where particles have an infinite inertia, the 
pressure term is thus dominant in the structure of the particle 
disk. 



<f>(z, u) = Ke 2tr -- e 



(17) 



where cr- is the initial thickness of the particle sub-disk and cr u 
is the initial (presumably small) velocity dispersion. The nor- 
malization constant K is related to the surface density X by 
L = 2nKcr z cr u . This initial condition assumes that the particles 
initially have a vanishing mean velocity, a velocity dispersion 
and are randomly distributed vertically, following a Gaussian 
distribution in this case. Other choices would lead to similar re- 
sults unless the symmetry with respect to the disk midplane is 
broken (particles initially located only in the upper half of the 
disk, for example). 

5. 1 . The corresponding moments 

Using definition (O, the density corresponding to the initial con- 
ditions dTvT > can be written as 



n(z, t) = Ke 2at exp 



( z 2 f 2 (t)) 






z 2 fl(t)\ 


{ 2<rf(r)J 


exp 







L(z,t) (18) 



where the exponential dependencies can be included into the 
Gaussian widths: 



4. General solution 

Equation (0 can be solved in the general case, using standard 
techniques (the method of characteristics for example). We will 
not develop this here as it is easier to check that the following 
solution is the general solution of Eq. 



ijf(z, u,f) = <p (zo(z, u, t), uq(z, u, f)) 



(13) 



where <f> is the initial condition and (zo, mq) is the initial phase 
space position of a particle located in (z, u) at a time t . The initial 
positions are given by 



zo=Mt)ze a < 
uo = Mt)z e a < 



f 2 (t)u e a < 
-f 3 (f)u e a >, 



(14) 



where f\ , and fa are time dependent functions given by 



f/i(O=ch(/3?)-fsh08O 



Mi) = gsh(00 



(15) 



Mi) = chOSf) + fsh(J3t) 
and a and (3 are parameters depending on the Stokes number S : 

(16) 



P = V(a+ l)(a- 1) = Va 2 - 1. 



Here, a is always real and positive, while /3 is either a positive 
real number (for small Stokes numbers), or a purely imaginary 
number (for large Stokes numbers). However, whatever ft, the 
time functions f\, /2, and are always real. 



(o- u = (T u e 
and L(z, t) is a function defined as 



(19) 



L(z. 
with 



•J —< 



( u 2 fHt)^ 



exp 



exp 



'2< 2 (0, 



, s fifiz hfoz 

p(z, t) = — ^-r + 



2u? 



2cr' 2 



e~ pu du (20) 



(21) 



Let us define 



A = 



(22) 



Then: 



f 

J -ex 



- pll du. 



(23) 



Here we keep the /^-dependence explicitly as it will be useful in 
the following. We have L(z, t) = L p if p takes the value given by 
equation d2Tl i. This expression can be now written as 



L p = F(p) + F(-p) 
where F(p) is the Laplace transform of a Gaussian: 



(24) 



F( P ) 



e 2^ e 



u du = A .^e^erfc . (25) 
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Therefore, using the identity erfc(x) + erfc(-x) = 2, we have 



, 4V 

L„ = A y2ne~ 



(26) 



Using the definition (O and a similar procedure, we get for 
the first-order moment: 



nu = Ke 2at exp 



( *fim 



exp 



H(z,t), 



where 



H(z,t) 



X+oo 
U 
OQ 



e~^e- pu du. 



(27) 



(28) 



Now Hp can be written as a function of F p , using basic proper- 
ties of the Laplace transform: 



H p = d p F(-p) - dpF(p) = -dpLp 



Hence, 



u = —- 



dpLp 



(29) 



(30) 



Applying the same procedure to definition (0), we get the 
following expression for the particle fluid pressure: 



d\L p I dpL p 
P - n n ' 



and using Eq. 
pressure: 



(31) 

-7> \ '-p 
we finally get the simple expression for the 



P(z, t) = n(z, t)A 2 (t). 



(32) 



Here A can be understood as the velocity dispersion of the 
particles, and its ex pression l22l generalizes the asymptotics of 
IGaraud et all d2004l) . 

It is interesting to note that the general pressure only depends 
on z through the density. The Euler equation for the particle fluid 
involves only gradients of the pressure so this gradient is in fact 
directly proportional to the density gradient. In other words, the 
particles seem to follow a barotropic equation of state. This be- 
havior see ms closely related to the assumption of constant stop- 
ping time (IGaraud et a l. 2004). 

5.2. Scaling and asymptotics 
5.2.1. Case /3 real (small particles) 

When /3 is real (and positive), a is large, it is easy to check that 
A is a decreasing function of time. Its initial value is 



^(0) = cr u , 
and the large time asymptotics are 



A(t -> +oo) 



P<r z 



-(a+j3)t 



(33) 



(34) 



Starting from a presumably very low value (the initial veloc- 
ity dispersion of particles), A decreases with time, asymptotically 
like an exponential. The location where the pressure decreases 
the least is the midplane, since the density tends to a Dirac func- 
tion when times goes to infinity, 



n(z = 0, t) 



KA V2^e 2ff ', 



and 



P(z = 0,t) = A 3 K^e 2a ' 



(35) 



(36) 



(37) 



Even in the midplane, the pressure quickly drops with time, 
starting fr om a low value. In this case, as has long been known 
(see e .g. ICuzzi et alJ Il993t iDubrulle et all ll995UGaraud et alJ 
2004), it is safe to neglect the pressure of the particle fluid. 
Consequently, particles steadily sediment towards the disk mid- 
plane at a velocity close to the terminal velocity W°° = OJ 2 ZT 
(defined as the asymptotic velocity of any body free-falling in a 
fluid), which is the asymptotic velocity of both the particle fluid 
and individual particles here. 

5.2.2. Case /3 imaginary (large particles) 

We define y such that /3 = iy. In such a case, A can be written as 

A = e~ M = (38) 

J (cos yt + 2 sin ytf + [^f ^sin 2 yt 

For small initial velocity dispersions, the location of the maxima 
of A corresponds to 

tanyf = -- (39) 

r 

and the amplitude of these maxima is 

KaAt) = a- J- Ja 2 +y 2 e~ a ' '. (40) 
a ' 

It is crucial to note here that the amplitude of the maxima of 
the time-dependent velocity dispersion A is independent of the 
initial velocity dispersion cr u . The velocity dispersion is created 
and maintained by the vertical gravity. 

Since the midplane density is still given by Eq. d35l l. we get 
the amplitude of the midplane pressure maxima: 



P max (t) = K^irZ yja 2 + y 2 . 



(41) 



These pressure maxima have a much lower decreasing rate than 
in the non-oscillating case, both because it involves only a here 
and because a has a lower value. 

The midplane pressure as a function of time for different 
values of the Stokes number S is displayed in Fig. Q] The 
transition between the case with oscillation S > 1/2 and the 
overdamped case S < 1/2 appears very clearly. The pressure is 
computed for a ratio 2* =0.1. This ratio defines the pressure at 
t — and the width of the maxima, but neither the amplitude of 
the maxima nor the decreasing rate. 

FigureQ]illustrates that, while the pressure of the solids fluid 
can be safely neglected for small values of the Stokes number 
S , the situation is drastically different when the Stokes number 
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Fig. 1. Pressure of the particle fluid in the midplane as a function 
of the dimensionless time (cjt), for different values of the Stokes 
number S 



reaches S =1/2. When particles oscillate around the midplane, 
the fluid macroscopic velocity is different from the terminal ve- 
locity vv 00 as the macroscopic pressure force slows down the ver- 
tical contraction of the particle fluid. In other words, the velocity 
of the particle fluid is not that of individual particles and the fluid 
has a nonvanishing internal energy (contrary to the case for small 
particles). 



6. Discussion and conclusions 

We have shown that, in general, the pressure of the particle fluid 
cannot be neglected as soon as their Stokes number is more than 
0.5. The particles seen as a fluid can have an internal energy, 
and the pressure force can partly balance gravity. In the extreme 
case where particles have an infinite inertia, the structure of the 
particle fluid is in fact hydrostatic. The existence of pressure is 
completely independent of the presence or absence of collisions 
between solid particles. Instead, the pressure is the result of a 
velocity dispersion. While velocity dispersion can be the result 
of collisions, it can also be created dynamically when external 
forces excite the velocity dispersion (gravity in the present case) 
or when the flow itself is converging locally. Even though our 
results concern only a very simple and specific problem, they 
have implications for a much broader variety of astrophysical 
and non-astrophysical flows. 

The Stokes number can be defined in a more general way as 
the ratio of a dynamical timescale, creating velocity differences 
between the solids and the underlying fluid (gas or liquid), and 
the stopping time, the characteristic timescale for the particles to 
reach the fluid velocity. There are many processes able to create 
velocity differences between the gas and solids, gravity being 
only one of them. For example, when two rivers join, particles 
embedded in the rivers do not see the confluence like the water 
does and this can create a velocity difference between water and 
the particles, and a velocity dispersion in the particle flow at the 
confluence. 



For particles embedded in a turbulent flow, it is common to 
define a scale-dependent Stokes number by S v = Tft v , where t n 
is the turnover time of a turbulent eddy of size rj. This Stokes 
number depends on the size of the turbulent eddy an d scales as 
S„ ex rj 1 ^ for Kolmogorov-like turbulence (see e.g. lCuzzi et al.l 
1200 ll) . The smallest reachable eddy is dependent on the strength 
of turbulence, quantified by the Reynolds number of the flow 
Re = —, where U and L are characteristic velocities and 
length scale, respectively, and v is the molecular kinematic 
viscosity of the fluid. For Kolmogorov turbulence, the smallest 
turbulent scale ^ mi „ of the flow scales as r],„i„ cc Re~ 3 ^ 4 . For 
flows with large enough Reynolds numbers (many astrophysical 
and geophysical flows reach Reynolds numbers over 10 10 ), the 
smallest turbulent scales have S n smaller than 1, even for small 
grains. This is a well-known behavior of particle-laden flows. 
Particles are centrifuged out of th e turbulent eddies when their 
scale Stokes number reaches unity (Squires & E atonll99ll) . This 
creates a converging flow of particles towards regions of low 
vorticity. Consequently, this induces a velocity dispersion inside 
the particle flow and, for similar reasons to those described in 
this paper, the particle fluid develops a pressure. As a side note, 
when turbulence modelling is considered, turbulent pressure 
terms can arise from fluc tuating or subgrid velocity dispersions 
dDobrovolskis et all 1 1999b . However, turbulent pressure terms 
are a consequence of Reynolds averaging or spatial filtering 
and disappear when turbulence is directly simulated instead of 
modelled. 

This is hopefully a good illustration of the care one has 
to take whenever a fluid equation is applied to a collection of 
solid particles in an astrophysical flow. Whenever any dynam- 
ical timescale becomes smaller than the characteristic coupling 
timescale between the particles and the underlying fluid, the par- 
ticles cannot be considered as a pressureless fluid, whatever the 
collision rate between particles. 
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